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0.  Introduction  and  Summary 

This  paper  is  concerned  with  methods  for  generating,  in  a  digital 
computer,  sequences  of  numbers  which  will  serve  as  realizations  of  a 
sequence  of  random  variables.  The  problem  has  two  parts  -  the  problem 
of  producing  independent  uniform  random  variables,  and  the  problem  of 
producing  arbitrary  random  variables  in  terms  of  uniform  random  variables . 
Section  1  describes  arithmetic  procedures  for  generating  uniform  variables. 
It  contains  nothing  new.  Section  2  considers  properties  of  linear  trans¬ 
formations,  modulo  1,  of  uniform  variables.  These  are  the  theoretical 
analogues  of  the  arithmetic  procedures  for  generating  uniform  variables. 

The  rest  of  the  paper,  sections  3-7,  is  devoted  to  the  problem  of  pro¬ 
ducing  arbitrary  random  variables  in  terms  of  uniform  random  variables. 

1 .  Producing  Uniform  Random  Variables. 

The  most  widely  used  current  method  for  generating  uniform  numbers 
is  to  successively  generate  the  residues  of  some  integer  by  a  linear 
transformation.  This  was  first  suggested  by  Lehmer,  [l] ,  and  has  be¬ 
come  the  standard  procedure.  As  a  simple  example,  in  the  ring  of 
residues  of  100,  the  linear  transformation  llx  +  7  starting  with,  . 
say  X  =  54,  generates  the  sequence 

54,1,18,5,62,89,86,53,... 

and  the  apparently  haphazard  succession  of  integers  sviggests  that  the 
sequence  might  well  serve  as  a  sequence  of  random  digits.  It  is  easy 


I 
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to  choose  the  modulus  and  the  linear  transformation  so  that  a  large 
number  of  residues  can  be  generated  without  repetition,  see  [2] -[6] 
for  details  and  further  references.  Rotenberg,  [5],  has  given  what 
is  probably  the  moat  suitable  method  for  digital  computers:  To  generate 
a  sequence  of  residues  put 

x^^^  =  (10^  +  l)x^  +  c,  modulo 

for  a  decimal  computer, with  k  ^  2  and  c  not  divisible  by  2  or  5,  and 
x^^^  =  (2*^  +  l)x^  +  c  modulo  2"' 

for  a  binary  computer,  where  k  ^  2  and  c  is  odd.  These  transforma¬ 
tions  generate  all  of  the  possible  residues  for  their  particular  modulus, 
and  are  particularly  suitable  for  digital  computers  since  no  proper 
multiplications,  only  shifts  and  additions,  are  required. 

As  an  example  of  how  the  procedure  works  in  a  computer,  suppose 
we  have  a  decimal  machine  which  manipulates  blocks  of  8  digits.  We 

g 

might  generate  the  residues  of  10  systematically  by  the  transformation 
x^^^  =  10^ x^  +  x^  +  231  (mod  10^) 

/ 

Starting  with,  say,  21329312,  we  would  have 

x^  =  21329312 
^  29312231 

X  =  50641543 
41543231 
X  =  92184774 

84774231 

X  =  76959005 
^  59005231 

X  =  35964236 


g 

Rather  than  viewing  the  x's  as  residues  of  10  , 


we  would 
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g 

think  of  them  as  fractions  of  10  ,  and  thus  the  above  sequence  would 
provide  us  with 

.21329312,  .50641543,  .92184774,  .76959005,  .35964336,... 

as  a  realization  of  the  sequence  u^,U2,u^,...  of  independent  uniform 
[0,1]  random  variables.  When  the  iterates  are  viewed  as  fractional 
parts  of  the  modulus ,  our  trans  formation  takes  the  form 

T(x)  =  lOOlx  +  .00000231  (modulo  l) 

and  this  leads  us  to  a  consideration  of  the  properties  of  such  linear 
transformations . 

2.  Linear  Transformations,  Modulo  1,  of  a  Uniform  Variable. 

Let  u  be  a  uniform  [0,l]  random  variable,  and  suppose  we 
produce  a  new  random  variable,  completely  dependent  on  u,  by  a 
linear  transformation  mod  1, 

T(u)  =  nu  +  a  (mod  1)  n  an  integer,  0  <  a  <  1. 

It  is  interesting  to  note  that  the  sequence 

(l)  u,T(u) ,T^(u)  ,.. . 

although  a  sequence  of  completely  dependent  variables,  still  has 
properties  similar  to  those  of  a  sequence  Uq,u^,U2,...  of  Independent 
uniform  [0,l]  random  variables. 

Suppose  we  look  at  the  joint  distribution  of  u,T(u) .  They  are 
individually  uniform  [0,l],  but  of  course  their  joint  distribution 
is  degenerate.  In  fact,  the  point  (u,T(u))  is  uniformly  distributed 
over  the  graph  of  T,  and  the  larger  n,  the  more  the  graph  of  T 
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tends  to  fill  the  unit  square, 
n  =  10  and  100  are  dravm. 


T(u)=10u(M0D  1) 


In  figure  1,  graphs  of  T  for 


0  1 

T(u)=100u(MOD  1) 


Thus,  for  large  n,  the  ped.r  u,T(u) ,  while  completely  dependent, 
might  well  serve  as  a  pair  of  independent  viniform  [0,l]  random  variables 

It  is  easy  to  show  that  the  sequence  (l)  satisfies  the  central 
limit  theorem;  in  fact  the  sum  u  +  T(u)  +•••+  T  (u)  quickly  converges 
to  the  distribution  of  the  analogous  sum  u^  +  u^  +•••+  Uj^  of 
independent  uniform  [0,l]  random  variables.  Figure  2  shows  the 
density  functions  of  u  +  T(u)  and  u  +  T(u)  +  T  (u)  compared  with 
those  of  Uq  +  u^  and  Uq  +  u^^  +  u^j  with  T(u)  =  lOu  (mod  l) . 


1.0 
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Density  function  of  u+T(u)  compared  with  that  of  Uq +01  . 


Density  function  of  u+T(u)+'I^  (u)  compared  with  that  of  Uo+Uj+Ua  . 

The  above  considerations  sxiggest  that  we  should  make  n  large  in 

2 

our  arithmetic  procedure  for  generating  u,T(u),T  (u) , . . . .  Coveyou,  [6], 

has  suggested  that  T  be  chosen  so  as  to  reduce  the  correlation  between 

u  and  T(u)  .  However,  if  n  is  large,  that  correlation  will  be 

negligible.  We  can  easily  find  the  correlation  between  u  and  T(u) . 

Suppose  we  use  an  overscore  to  indicate  the  fractional  part  of  a 

number;  thus  T(u)  =  a  +  nu  and  in  general,  x  =  x  -  [x] .  Now  we  may 
d 

write  u  =  — +  —  where  d,  is  a  random  integer  taking  values  0,l,...,n-l 
n  n  1 

with  probabilities  and  v  is  \iniform  [0,l],  independent  of  d^. 

Then  T(u)  =  a  +  nu  =  a  +  d^  +  v  =  a  +  v.  It  is  easy  to  verify  that 
a  +  V  is  uniform  [0,l]  and  that  the  expected  value  of  v(a  +  v)  is 
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Thus 


nE[uT(u) ]  =  E[(d^+v) (a+v) ]  =  E(d^)E(a+v)  +  E[va+v], 
and  it  easily  follows  that  the  correlation,  p,  between  u  and  T(u)  is 


P  = 


6(a  -  ^)  -  5 


n 


Thus  “  ^  <  P  p  =  0  when  a  =  ^  +  ^'s/3*  But  in 


2n  '  '  n 

p  will  be  small  when  n  is  laiige. 


any  case, 


3*  Producing  Arbitrary  Random  Variables 

If  we  are  willing  to  grant  the  adequacy  of  arithmetic  procedures 
for  generating  uniform  numbers,  then  we  may  assume  we  have  a  sequence 

of  independent  [0,l]  random  variables  and  turn  to  the  problem  of 
producing  arbitrary  random  variables  in  terms  of  the  u's.  An  obvious 
way  of  generating  a  random  variable  x  with  continuous  distribution 
F  is  to  put  X  =  F~^(u) ,  but  this  is  not  usually  feasible  in  a 
computer  -  it  may  take  too  long  to  evaluate  F~^(u)  in  a  subroutine, 
while  a  direct  table  look-up  may  require  too  many  storage  locations 
for  the  necessary  precision.  The  storage  requirement  may  be  reduced 
by  using  a  coarser  mesh  and  interpolation,  but  even  linear  interpola¬ 
tion  usually  requires  one  multiplication,  and  we  are  concerned  with 
procedures  in  which  the  time  for  a  single  multiplication  would  be  a 
preponderant  part  of  the  average  running  time. 


A  general  procedure  which  will  provide  fast  programs  for  generating 
X  with  distribution  F  may  be  based  on  representing  F  as  a  mixture. 
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The  method  Is  described  In  [?]>  from  which  this  s\nnmary  Is  taken. 

The  Idea  Is  roughly  as  follows t  Suppose  we  have  a  method  for 

providing  a  random  variable  with  distribution  function  and 

method  takes  10  units  of  time.  Suppose  we  also  have  a  method 

M2  for  providing  a  random  variable  with  distribution  Ainctlon 

G2>  and  method  M2  takes  500  units.  If  we  can  represent  F  as  a 
mixture  of  G^  and  G2>  say 

F(a)  =  .99G^(a)  +  .01G2(a), 

then  we  may  use  Uj^  ,U2 lU^ » .  •  •  to  provide  a  random  variable  x  with 
distribution  F  as  follows:  If  u^  <  .99>  we  use  method  M^  on 
U2,U2,...  to  generate  y^^  and  put  x  =  y^^.  If  .99  < 
use  method  M2  on  to  generate  y2  and  put  x  =  y2.  The 

average  running  time  will  then  be  (.99)10  +  (.01)500  =  14.9  units,  plus 
the  time  necessaiy  to  test  u^^  <  .99. 

This  illustrates  the  basic  principle  of  the  simple  device  which 
provides  programs  with  very  short  average  running  times  -  we  represent 
F  as  a  mixture  of  distributions, 

(1)  .  F(a)  =  p^G^(a)  +  p2G2(a) , 

in  such  a  way  that  p^  is  close  to  1  £uid  the  time  to  generate  a 
random  variable  y^  with  distribution  G^  is  small;  then  most  of 
the  time  we  put  x  =  y^^.  Even  though  G2,  the  correcting  distribution, 
may  be  quite  complicated  and  difficult  to  handle,  we  still  have  a 
short  average  running  time,  since  G2  must  be  handled  so  infrequently. 

In  searching  for  representations  of  F  such  as  (l) ,  if  we  have 
a  G^  which  shows  promise,  then  we  find  the  largest  p^^  so  that 
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f(x)  -  p^gj^(x)  is  non-negative.  Furthermore,  If  we  make  g^^  a  . 
rectangle  or  a  mixture  of  rectangles,  we  can  expect  to  have  very 
short  running  times,  especially  if  the  rectangles  are  chosen  so  as 
to  exploit  the  particular  features  of  the  computer  in  question. 

References  [7], [8]  give  details  of  how  this  method  may  be  applied.' 
The  general  idea  is  to  break  the  density  ftmction  into  two  parts, 
rectangles  auid  "teeth”,  as  for  example.  In  this  figure: 


Then  a  random  variable  with  density  g^  may  be  produced  quite  rapidly; 
only  occasionally  will  the  computer  have  to  deal  with  the  density  gg, 
and  even  then,  a  reasonable  procedure  may  be  constructed  along  these 
lines:  One  of  the  teeth  of  is  chosen  with  appropriate  probability, 

then  a  random  variable  with  the  nearly-linear  density  function  represented 
by  that  tooth  must  be  generated.  We  will  describe  a  modified  rejection 
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technique  for  generating  such  a  variable  In  the  next  section. 

4.  Generating  a  Random  Variable  With  a  Nearly  Linear  Density. 

Suppose  we  want  to  generate  a  random  variable  z  with  a  nearly 
linear  density  function  g(x),  0  <  x  <  c,  such  as  In  this  figure: 


Let  AD  euid  flC  be  parallel  segments  enclosing  g(x).  Then  z  may 
be  produced  as  follows: 

1)  Choose  independent  uniform  [0,l]  random  variables  u  and  v, 
let  M  =  max(u,v)  ,  m  =  min(u,v)  . 

2)  Choose  a  point  (x,y)  uniformly  from  triangle  OBC,  for 
example,  by  putting’^ 

X  =  cm 
y  =  bM  -  bm 

3)  If  y  ^  -  b^^^c  +  a,  i.e.,  if  (x,y)  is  in  triangle  QAD, 

put  z  =  X  and  stop. 

4)  If  -  hx/c  +  a  <  y  <  g(x),  put  z  =  x  and  stop. 

5)  If  g(x)  <  y,  return  to  step  1  and  try  a  new  pair  u,v. 

*In  fact,  any  one  of  the  six  representations  (c-cM,bm)  ,(cm,b-bM)  , 

(cM-cm,bm)  ,(cm,bM-bm)  ,(c-cM,bM-bm)  ,(cM-cm,b-bm)  will  provide  a  point 
(x,y)  miformly  distributed  over  triangle  OBC.  See,  e.g.,  [9]. 


If  we  can  enclose  g(x)  in  a  narrow  enough  atrip,  we  will  produce 
X  most  of  the  time  without  computing  g(x)  -  step  3  will  provide  x 
with  probability  near  1.  Note  that  the  test  in  3) ,  y  ;^-bx/c  +  a,  Is 
equivalent  to  bM-bm  ^  -  bm  +  a,  or  M  ^  a/b,  and  a/b  should  bo 
close  to  1  if  g(x)  is  nearly  linear.  We  may  summarize  the  pro¬ 
cedure  as  follows: 


renerate  a  random  variable 


with  a  nearly  linear  densitv 


function  g(x) ,  O^x^c, 


such  as  this  or  this 


like  this 


or  this 


Then: 

1.  Choose  independent  uniform  [0,l]  random  variables  u  and  v. 

2.  ^  max(u,v)  <  a/b  put  z  =  c  inin(u,v)  and  stop. 
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3.  If  not,  test;  b  |u  -  v|  <g[c  inin(u,v)]?  If  yea .  put 
z  =  c  mln(u,v)  j  ^  iw,  ^  step  1  and  try  again. 

Note  that  for  g  concave  up,  as  in  the  left  case  above,  we  may  put 

b  =  g(0)  and  a  =  bx^y^c  +  gCx^),  where  g'Cx^)  =  -  b/c.  If  g  is 
concave  down,  such  as  in  the  rlg;hib  case,  then  a  =  g(0)  and 
b  =  -  eg '(c)  will  do. 

5.--  Special  Methods  for  Particular  fiandom  Variaibles. 

It  is  sometimes  possible  to  take  advantsge  ef  ^Mitdcular  properties 
of  the  desired  random  variable  in  order  ito  .geaerat®  it.  For  example, 
as  mentioned  in  [lO] ,  an  expaaeiatial  ramdmm  -warlable  X,  density 
e  ,  X  >  0,  may  be  generated  by  gmaitltiiiig 

X  =  m  +  min(u^,U2, .  • .  ,u^) 

where  m  and  n  are  random  integers  with  distributions: 

P[m  =  k]  =  (e  -  l)e’^'*'^  k  =  0,1,2,... 

P[n  —  k]  —  -~iy  ^  ~  1,2,3,.». 

A  possible  method  for  generating  a  number  of  exponential  random 

variables  simultaneously  is  as  follows:  Let  n  be  fixed,  let 

n  =  (Pi>P2»  •  •  *  »Pn+l^  viniformly  distributed  over  the  regular  simplex 

-X  n 

®BriH  ’  * ’"’’^n+l"^^  ^  density  ®-^j-,x  >  0, 

and  be  independent  of  tt.  Then  the  n+1  random  variables  zp^,zp2»  • .  •  ,zp^^j^ 
are  mutually  independent  exponential  random  variables.  A  proof  is  in  [9]. 

We  may  produce  (p^,P2,  • . .  uniformly  over  by  ordering  n 

independent  uniform  [0,l]  random  variables  ^^(2)  ^***— ^(n)* 
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then  putting  =  '*(l)*P2  "  ^^(2)  "  “(l)*P3  “  “(3)  "  ^(2)»***»Pn 
=  '"(n)  -  '"(n-1)  'Pn+1  =  ^  "  "(n)  * 

A  method  for  generating  a  pair  (x,y)  of  normal  random  variables 
may  be  based  on  the  polar  representation,  putting 


X  =  p  cos  9 


y  =  p  sin  0 

r  2  2 

where  0  Is  uniform  [0,2tt]  and  p  Is  distributed  as  vx  +  y  ,  i.e., 

2  2/ 

^  has  the  exponential  distribution:  P[p  <  a]  =  1  -  e  '  .  As  pointed 

out  by  Von  Neumann,  [ll] ,  cos  0  and  sin  0  may  be  generated  directly  in 

the  form  ^2'  *  ^2^^1  ^  ^2'  where  and  are  uniform 

2  2 

[-1,1],  conditioned  by  ^  ^2  —  alternatively,  avoiding  the 


square  root,  in  the  form 


^^1^2 
2^2 
^1  +  ^2 


2  2 
^1  -  ^2 
2  2 
^1  +  ^2 


This  procedure  is  quite  slow  for  general  use,  but  it  does  provide  a 
convenient  way  of  handling  the  tail  of  the  normal  distribution. 


6.  Generating  a  Discrete  Random  Variable. 

In  section  7  we  will  describe  a  general  procedure  for  generating 
arbitrary  random  variables.  It  is  based  on  the  ability  to  generate  a 
discrete  random  variable  in  a  quick  and  simple  manner.  If  x  is  a 
discrete  random  variable  with  P[x  =  a^]  =  p^,  then  an  obvious  method 
for  generating  x  in  a  computer  is  to  generate  a  uniform  [0,l]  random 
variable  u  and  put  x  =  a^  if 


+•••+  Pj^_^  <  u  <  Pj^  +•••+  p^. 


Various  search  techniques  based  on  this  method,  (see,  e.g.,  [12]),  may 
be  used,  but  they  often  are  relatively  complicated  progran^  that  take 
too  long.  We  will  outline  a  procedure  that  leads  to  short  and  fast 
programs.  To  avoid  notational  difficulties,  consider  this  particular 
example : 


value  of  X  probability 


a 

b 

c 

d 

e 

f 

g 

h 

i 

3 

k 

m 


.023  =  .Df.02+.003 
.038  =  .0f.03+.008 
.074  =  .0t.07+.004 
.103  =  .l+.00f.003 
.148  =  .1+.04+.008 
.206  =  .2+.0Ot.006 
.140  =  .1+.04+.000 
.101  =  .l+.OOf.OQl 
.093  =  .0f.09+.003 
.037  =  .0f.03+.007 
.026  =  .0f.02+.006 
.011  =  .Of.Ol+.OOl 
.35  .050 


We  have  represented  the  probabilities  in  decimal  form  and  will 
consider  a  decimal  computer  -  the  ideas  apply  equally  well  to  a  binary 
machine.  Suppose  our  decimal  computer  has  storage  blocks  (words) 
which  may  be  called  for  by  number,  as  most  do.  Then  the  fastest  method 
for  generating  the  discrete  random  variable  x  is  probably  this: 


In  memory  locations  0-999,  store  23  a's,38  b's,74  c's,103  d's, 
•••,  11  m's.  Then  if  u  =  .dj^d2d2d^. . .  is  a  uniform  random  variable, 
generated  in  the  computer  in  some  way,  look  up  the  number  in  location 
d^d2d2  and  let  that  be  x. 


In  effect,  we  have  an  um  containing  1000  balls  -  23  marked  a, 

38  marked  b,  •••  euid  we  choose  one  of  the  1000  at  random.  While  this 


is  the  fastest  method,  it  requires  1000  storage  locations.  We  will 
consider  an  improvement  on  the  fastest  method.  It  is  an  improvement 
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In  that  it  uaeB  much  less  memory  space  and  takes  little  longer^  To 
continue  the  um  analogy,  suppose  we  fill  3  urns  according  to  the 
digits  of  the  probabilities  for  x  as  follows: 

um  1-1  d,l  e,2  f's,l  g,l  h  (first  digit) 

um  2-2  a's,3  b's,?  c’s,...,l  m  (second  digit) 

um  3-3  a's,8  b's, 4  c's,3  d's,***,l  m  (third  digit) 

Now  choose  um  1  with  probability  .600,  um  2  with  probability  .350, 
um  3  with  probability  .050)  then  choose  a  ball  at  random  from  that 
um.  We  then  require  only  91  balls  instead  of  1000,  and, on  the  average, 
only  a  few  more  steps  than  the  fastest  um  scheme. 


In  a  computer,  we  would  "stack"  the  urns  -  um  1  would  occupy 
storage  locations  0-5,  um  2  locations  6-40,  and  um  3  locations 
41-90.  An  outline  of  a  program  for  generating  this  particular  x 
might  run: 


Load  constants  in  memory  locations 


Lo^tion  Contents 


0  d-* 

10  b 

20  e 

30 

i 

40  m 

1  e 

11  c 

21  e 

31 

i 

41  a 

2  f 

12  c 

22  g 

32 

1 

42  a 

3  f 

13  c 

23  g 

33 

i 

43  a 

4  g 

14  c 

24  g 

34 

i 

U  b 

5  h 

15  c 

25  g 

35 

i 

45  b 

6  a 

16  c 

26  i 

36 

i 

46  b 

7  a 

17  c 

27  1 

37 

i 

47  b 

8  b 

18  e 

28  1 

38 

k 

48  b 

9  b 

19  e 

29  1 

39 

k 

49  b 

Let 

u  =  .dj^d 

2*^3*  •* 

be  a 

uniform 

0-90  according  to  this  scheme: 


50 

b 

60 

e 

70 

f 

80 

i 

90  m 

51 

b 

61 

e 

71 

f 

81 

i 

91 

52 

c 

62 

e 

72 

f 

82 

i 

92 

53 

c 

63 

e 

73 

h 

83 

•j  ■ 

93' 

54 

c 

64 

e 

74 

i 

84 

k 

94 

55 

c 

65 

e 

75 

i 

85 

k 

95 

56 

d 

66 

e 

76 

i 

86 

k 

96 

57 

d 

67 

f 

77 

i 

87 

k 

97 

58 

d 

68 

f 

78 

i 

88 

k 

98 

59 

e 

69 

f 

79 

J 

89 

k 

99 

1]  random  variable.  Let  a(n) 
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be  the  contents  of  memory  location  n.  Then 


1) 

If 

d^  <  6 

put 

X  = 

a(dj^) 

2) 

If 

60  ^  d^d2  <  95 

put 

X  = 

a(did2  - 

54) 

3) 

If 

950  1  d^d2d3 

put 

X  = 

a(did2d3 

-  919) 

Examples 

• 

• 

u  =  .217... 

X  =  a 

(2)  = 

=  f 

u  =  .728... 

X  =  a 

(18) 

=  e 

u  =  .963.  •  • 

X  =  a 

(U) 

=  b 

7.  The  Latest  and  Best  Method 

The  ability  to  provide  a  discrete  random  variable  quickly  and 
easily  serves  as  the  basis  for  a  general  method  for  producing  arbitrary 
random  variables.  This  general  method  is  very  fast  and  easy  to  program 
for  a  computer.  It  is  the  latest  and  beat  of  a  series  of  procedures 
we  have  recently  considered.  The  idea  is  quite  simple  -  to  produce  a 
random  variable  y  with  distribution  F,  we  write  y  =  x  +  v,  where 
X  is  discrete  and  v  uniform.  This  amounts  to  replacing  the  distri¬ 
bution  F  by  a  piecewise-linear  distribution  G.  By  making  x  a 
random  multiple  of  a  sma].l  increment  6,  we  can  make  G  as  close  to 
F  as  we  please.  Details  of  the  procedure  for  the  common  distributions 
will  appear  elsewhere.  We  illustrate  the  procedure  by  showing  how  to 
generate  an  exponential  variable:  To  generate  y  with  P[y  <  a] 

=  F(a)  =  1  -  e”®^,  we  generate  y  in  the  Interval  0  <  y  ^  4*6  by 
putting  y  =  X  +  V  where  x  takes  values  0, .1, .2, . . . ,4.5  with 

_  Jl  _  iiti 

P[x  =  =pj^=e  -e  k  =  0,1,2, ...  ,45, 

and  V  is  uniform  [0,.l].  Values  of  y  in  the  Interval  4*6  ^  y  <  “ 
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are  produced  by  taking  the  logarithm  of  a  vmlform  number. 

The  probabilities  for  x  are  expressed  only  to  enough  decimal 
places  to  give  the  accuracy  necessary  for  the  portion  of  F  to  vhlch 
Pj^  applies,  and  the  p’s  are  rounded  off,  high  or  low,  In  such  a  way 
that  the  cumulative  distribution  p^  +  Pg  +•••+  pj^  stays  close  to 
F.  The  p’s  are  displayed  In  the  table  of  constants  for  loading  the 
memory;  the  digits  of  the  p's  are  used  to  provide  the  compact  storage 
method  described  In  section  6. 

The  procedure  below  will  produce  a  value  y  with  distribution  G 
for  which  |f~^  -  G~^ |  <  .001,  that  Is,  It  will  return  a  value  y 
which  differs  from  the  "true"  value  by  less  than  .001.  The  steps  of 
the  procedure,  very  easy  to  program  for  a  computer,  are  as  follows: 

To  produce  an  exponential  random  variable  Y,  density  e~^, 
y  >  0,  store  the  constant  a(n)  ,  n  =  0,1, . . .  ,569  iin  memory  location 
n.  Let  u  =  .d^dgd^. . .  ^  a  imlform  [0,l]  random  variable ,  the 

d's  its  decimal  digits.  Then 


1) 

If 

0 

^did^ 

<  79 

put 

Y  —  a(d^d2)  t  *  *  * 

2) 

If 

790 

<  d^d2d3 

<  969 

put 

Y  =  a(d^d2d3  -  711)  +  -Od^d^dg  ... 

3) 

If 

9690 

1  d^d2d3d^ 

<  9888 

put 

Y  =  a(d^d2d3d^  -  9432)  +  .Oigd^dg  .. 

4) 

If 

98880  <  d,d_d_d,dc 

<  98994 

put 

Y  =  a(d^d2d3d^d5  -  98424)  +  -Odgd^dg 

5) 

If 

98994  1  dj^d^d^d^d^ 

put 

Y  =  -  ln(u  -  .98994) 

The  constants  for  loading  the  memory  are  given  in  the  table. 

The  digits  of  the  probabilities  for  x  provide  the  frequencies  of  the 
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values  to  be  stored.  The  second  digit  provides  the  frequencies  for 
table  1,  the  third  digit  the  frequencies  for  table  2,  etc.  The  tables 
are  then  "stacked"  in  the  memory. 


Value 

.0 

.1 

.2 

.3 

.4 

.5 

.6 

.7 

.8 

.9 

1.0 

1.1 

1.2 

1.3 

1.4 

1.5 

1.6 

1.7 

1.8 

1.9 
2.0 
2.1 
2.2 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 
3.0 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6 

3.7 

3.8 

3.9 
4.0 

4.1 

4.2 

4.3 

4.4 

4.5 


Table 
12  3  4 


Constants  for  Loading  the  Memory  in  Generating 
Exponential  Random  Variables 


.0  9 
.0  8 
.0  7 
.0  7 
.0  6 
.0-5 
.0  5 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 
.0 


Table 


Contents 


1  nine  O’s, eight  .l*s, seven  .2's,...,one  2.2 

2  five  0*s,six  .I's, eight  .2's,...,one  4*5 

3  two  0*s,four  .3*3, nine  .4*s,...,one  4.4 

4  eight  2.3’s,four  2.4's,eight  2.6's , . . . ,four  4.5's 

Table  Memory  Locations 

1  0-78 

2  79-257 

3  258-455 

4  456-569 


Size  of 
Table 

79 
179 
198 
114 


8 

4 
0 
8 
8 
0 

3 

5 
7 
9 
7 
2 

6 

1 

4 

4 
2 

5 

7 
4 

8 
8 


.00104 
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In  conclusion,  It  appears  that,  except  for  certain  rare  situations 
In  which  special  tricks  such  as  those  In  sections  3  or  5  may  be  more 
suitable,  the  best  method  for  generating  a  random  variable  y  Is  to 
generate  a  oonlform  niimber  u  =  .d^^d^d^...  by  one  of  the  arithmetic 
procedures,  with  enough  digits  In  u  that  the  first  few  may  be  used 
to  assign  a  value  to  the  discrete  variable  x  by  the  compact  storage 
method  of  section  6,  then  put  y  =  x  +  v,  where  v  is  formed  from  the 
digits  in  the  latter  part  of  u.  Occasionally,  if  u  is  close  to  1, 
the  tall  of  the  distribution  will  provide  y,  rather  than  the  simple 
procedure  y  =  x  +  v.  The  compact  storage  procedure  for  providing  x 
makes  the  method  very  easy  to  program,  and  the  average  time  to  produce 
y  is  little  longer  than  the  time  to  make  a  comparison,  look  up  a  stored 
value,  and  adjoin  it  to  the  last  digits  of  u.  Note  that  forming  x  +  v 
is  not  really  an  addition,  but  merely  a  juxtaposition  of  the  digits  of 
x  and  those  of  v. 

The  average  time  to  produce  u  in  one  of  the  common  U.S.  computers, 
the  IBM  7090,  is  about  30  x  10~^  seconds ,  and  the  average  time  to  produce 
y  from  u  is  about  the  same,  so  that  random  variables  yi>y2’**’ 
be  produced  at  the  rate  of  about  17000  per  second. 

The  status  of  the  problem  of  generating  random  variables  in  a  computer 
is  then  this:  It  is  now  fairly  easy  to  develop  a  scheme  for  storing 
constants  in  a  computer  in  such  a  way  that  an  arbitrary  random  variable 
can.be  produced  quickly  and  easily.  The  instructions  for  the  computer 
are  easy  to  fomailate,  and  most  of  the  time  only  basic  operations  such 
as  testing  magnitudes,  shifting,  and  calling  for  stored  values  are 


required.  The  procedures  are  easy  enough  that  each  computer  center 
can  develop  standard  subroutines  of  their  ovm  for  providing  the  common 
random  variables.  Thus,  mathematicians  interested  in  using  artificial 
sequences  of  random  variables  to  guide  their  research  activities  can 
turn  their  attention  to  the  problem  of  what  to  do  with  the  large 
numbers  of  random  variables  that  the  computer  can  quickly  and  easily 
produce.  Their  first  consideration  will  of  course  be  to  judge  the 
suitability  of  the  sequence  for  their  particular  need.  In  this  connec¬ 
tion,  the  traditional  statistical  tests  on  the  random  digits  of  the 
\iniform  nimbers  are  of  limited  value.  More  appropriate  is  to  devise 
a  test  problem  similar  to  the  one  in  question  and  try  the  sequence  for 
that  problem  -  if  a  problem  on  the  random  packing  of  spheres,  try  the 
sequence  on  a  problem  on  the  random  packing  of  cubes  for  which  the  answer 
can  be  found.  If  one  wants  the  distribution  of  a  certain  statistic  for 
large  values  of  the  sample  size  n,  try  the  artificial  sequence  on  the 
largest  n  for  which  the  distribution  can  be  found  by  analytical  means. 

If  one  wants  to  find  the  distribution  of  the  length  of  a  chain  of  molecules 
joined  according  to  some  chance  mechanism,  test  the  artificial  sequence 
by  generating  a  related  random  walk  which  has  a  known  solution.  If  the 
problem  involves  some  difficult  function  of  order  statistics,  see  if  the 
artificial  sequence  gives  results  consistent  with  the  distributions  of 
some  of  the  explored  functions  of  order  statistics,  etc. 
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